Kevin Mader
6 March 2014
to a single, discrete value (usually true or false, but for images with phases it would be each phase, e.g. bone, air, cellular tissue)
2560 x 2560 x 2160 x 32 bit = 56GB / sample \[ \downarrow \]
2560 x 2560 x 2160 x 1 bit = 1.75GB / sample
Start out with a simple image of a cross with added noise
\[ I(x,y) = f(x,y) \]
The intensity can be described with a probability density function
\[ P_f(x,y) \]
By examining the image and probability distribution function, we can deduce that the underyling model is a whitish phase that makes up the cross and the darkish background
Applying the threshold is a deceptively simple operation
\[ I(x,y) =
\begin{cases}
1, & f(x,y)\geq0.5 \\
0, & f(x,y)<0.5
\end{cases} \]
While scalar images are easiest, it is possible for any type of image
\[ I(x,y) = \vec{f}(x,y) \]
The intensity can be described with two seperate or a single joint probability density function
\[ P_{\vec{f}\cdot \vec{i}}(x,y), P_{\vec{f}\cdot \vec{j}}(x,y) \]
A threshold is now more difficult to apply since there are now two distinct variables to deal with. The standard approach can be applied to both
\[ I(x,y) =
\begin{cases}
1, & \vec{f}_x(x,y) \geq0.25 \text{ and}\\
& \vec{f}_y(x,y) \geq0.25 \\
0, & \text{otherwise}
\end{cases} \]
This can also be shown on the joint probability distribution as
Given the presence of two variables; however, more advanced approaches can also be investigated. For example we can keep only components parallel to the x axis by using the dot product.
\[ I(x,y) =
\begin{cases}
1, & |\vec{f}(x,y)\cdot \vec{i}| = 1 \\
0, & \text{otherwise}
\end{cases} \]
We can tune the angular acceptance by using the fact \[ \vec{x}\cdot\vec{y}=|\vec{x}| |\vec{y}| \cos(\theta_{x\rightarrow y}) \] \[ I(x,y) = \begin{cases} 1, & \cos^{-1}(\vec{f}(x,y)\cdot \vec{i}) \leq \theta^{\circ} \\ 0, & \text{otherwise} \end{cases} \]
Going beyond vector the possibilities for images are limitless. The following shows a tensor plot with the tensor represented as an ellipse.
\[ I(x,y) = \hat{f}(x,y) \]
Once the variable count is above 2, individual density functions and a series of cross plots are easier to interpret than some multidimensional density hypervolume.
Ideally we would derive 3 values for the thresholds based on a model for the composition of each phase and how much it absorbs, but that is not always possible or practical.
For this exercise we choose arbitrarily 3 ranges for the different phases and perform visual inspection
The relation can explicitly be written out as
\[ I(x) =
\begin{cases}
\text{Void}, & 0 \leq x \leq 0.2 \\
\text{Clay}, & 0.2 < x \leq 0.5 \\
\text{Rock}, & 0.5 < x
\end{cases} \]
The implementations of basic thresholds and segmentations is very easy since it is a unary operation of a single image \[ f(I(\vec{x})) \] In mathematical terms this is called a map and since it does not require information from neighboring voxels or images it can be calculated for each point independently (parallel). Filters on the other hand almost always depend on neighboring voxels and thus the calculations are not as easy to seperate.
The simplist is a single threshold in Matlab:
thresh_img = gray_img > thresh
A more complicated threshold:
thresh_img = (gray_img > thresh_a) & (gray_img > thresh_b)
thresh_img = gray_img > thresh
thresh_img = map(lambda gray_val: gray_val>thresh,
gray_img)
boolean[] thresh_img = new boolean[x_size*y_size*z_size];
for(int x=x_min;x<x_max;x++)
for(int y=y_min;y<y_max;y++)
for(int z=z_min;z<z_max;z++) {
int offset=(z*y_size+y)*x_size+x;
thresh_img[offset]=gray_img[offset]>thresh;
}
bool* pix = malloc(x_size*y_size*z_size * sizeof (bool));
for(int x=x_min;x<x_max;x++)
for(int y=y_min;y<y_max;y++)
for(int z=z_min;z<z_max;z++) {
int offset=(z*y_size+y)*x_size+x;
thresh_img[offset]=gray_img[offset]>thresh;
}
| Radius | Mean Intensity | Sd Intensity |
|---|---|---|
| 2.000 | 0.0311 | 0.1623 |
| 2.938 | 0.0677 | 0.2370 |
| 3.875 | 0.1177 | 0.3061 |
| 4.812 | 0.1822 | 0.3688 |
| 5.750 | 0.2601 | 0.4196 |
| 6.688 | 0.3511 | 0.4567 |
| 7.625 | 0.4569 | 0.4763 |
| 8.562 | 0.5761 | 0.4690 |
| 9.500 | 0.7073 | 0.4259 |
Returning to the original image of a cross
We can now utilize information from neighborhood voxels to improve the results. These steps are called morphological operations.
Like filtering the assumption behind morphological operations are
Therefore these imaging problems can be alleviated by adjusting the balance between local and neighborhood values.
A neighborhood consists of the pixels or voxels which are of sufficient proximity to a given point. There are a number of possible definitions which largely affect the result when it is invoked.
The neighborhood is important for a large number of image and other (communication, mapping, networking) processing operations:
For standard image operations there are two definitions of neighborhood. The 4 and 8 adjacent neighbors shown below. Given the blue pixel in the center the red are the 4-adjacent and the red and green make up the 8 adjacent.
Formally these have two effectively equivalent, but different definitions.
In 3D there is now an additional group to consider because of the extra-dimension
If any of the voxels in the neighborhood are 0/false than the voxel will be set to 0
If any of the voxels in the neigbhorhood are 1/true then the voxel will be set to 1
Taking an image of the cross at a too-high threshold, we show how dilation can be used to recover some of the missing pixels
Taking an image of the cross at a too-low threshold, we show how erosion can be used to remove some of the extra pixels
An erosion followed by a dilation operation
A dilation followed by an erosion operation
Taking an image of the cross at a too-low threshold, we show how opening can be used to remove some of the extra pixels
Taking an image of the cross at a too-high threshold, we show how closing can be used to recover some of the missing pixels
A segmented slice taken from a cortical bone sample. The dark is the calcified bone tissue and the white are holes in the image
Applying the closing with even larger windows will close everything but being to destroy the underlying structure of the mask
We can characterize this by examining the dependency of the closing size and the volume fraction (note scale)
Using a better kernel shape (circular, diamond shaped, etc) the artifacts from morphological operations can be reduced
We can characterize this by examining the dependency of the closing size and the volume fraction (note scale)
We see that the dilation and erosion affects are strongly related to the surface area of an object: the more surface area the larger the affect of a single dilation or erosion step.